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Abstract 

Explicit second-order accurate finite-difference schemes for the 
approximation of hyperbolic conservation laws are presented. These schemes 
are nonlinear even for the constant coefficient case. They are based on 
first-order accurate upwind schemes. Their accuracy is enhanced by locally 
replacing the first-order one-sided differences with either second-order one- 
sided differences or central differences or a blend thereof. The appropriate 
local difference stencils are selected such that they give TVD schemes of 
uniform second-order accuracy in the scalar, or linear systems, case. Like 
conventional TVD schemes, the new schemes avoid a Gibbs phenomenon at 
discontinuities of the solution, but they do not switch back to first-order 
accuracy, in the sense of truncation error, at extrema of the solution. The 
performance of the new schemes is demonstrated in several numerical tests. 



INTRODUCTION 


The notion of TVD goes back to Harten's ground breaking paper [2]. He 

set out to construct second-order accurate finite-difference schemes for 

hyperbolic conservation laws which do not exhibit the spurious oscillations 
near discontinuities of the solution, as generated by the more classical 
second-order methods, such as Lax-Wendroff [7]. Harten discovered that any 
TVD scheme is a monotonicity preserving scheme which always produces 

oscillation-free weak solutions of hyperbolic problems. As Van Leer did in 
his early works [14,15], Harten realized that monotonicity preserving schemes 
can be only extended to higher-order accuracy when they are nonlinear even for 
the constant coefficient case. Following Roe [9] and Sweby [12], such 

nonlinear schemes can be interpreted as first-order accurate upwind schemes 
whose excessive numerical dissipation is counterbalanced by adding just enough 
antidiffusion to gain resolution while still suppressing pre- and post-shock 
oscillations. The amount of antidiffusion is monitored by flux limiters. 
Yet, at extrema of the solution, all formulations suggested so far for these 
limiters [1,9,12,14,15] restore the baseline first-order scheme by setting the 
antidif fusive flux to zero. This perpetual damping of extrema reduces the 
global error of TVD schemes to be of only first order in the L norm [4], 

oo 

which impedes particularly the approximate solutions for transient wave- 
tracking problems . 

To remedy this shortcoming, Harten and Osher [4] recently introduced 
their '’Uniformly high-order accurate NOnoscillatory (UNO) schemes." When 
applied to nonoscillatory initial conditions, UNO schemes preserve their 
nonoscillatory nature although they are allowed to accentuate local extrema- 



Thus, they are not only no longer TVD, but they also do no longer ensure that 
the numerical solution converges to the physically correct one [2,4,12]. The 
question whether one really has to sacrifice the TVD property for the sake of 
uniform accuracy motivated the present study. 

This paper describes explicit high resolution schemes of uniform second- 
order accuracy which are also TVD in the scalar case. The next section 
briefly recapitulates the Riemann problem and Roe's [8] (first-order accurate) 
approximate Riemann solver to set the stage for the construction of those high 
resolution schemes. As usually done, first the rather watered-down special 
case of scalar conservation laws is examined which gives us a class of new TVD 
schemes of Uniform second-order Accuracy (in short: TVDUA schemes). The 
design principles of those TVDUA schemes are then applied to hyperbolic 
systems using a field-by-field decomposition. Eventually, the TVDUA schemes 
and their high-resolution counterparts for hyperbolic systems (which are not 
necessarily TVD for nonlinear systems) compete with conventional TVD schemes 
in several comparative numerical studies, where they show their superior 
accuracy at virtually no computational extra cost. 


THE PROBLEM AND ITS FIRST-ORDER ACCURATE APPROXIMATE SOLUTIONS 


Consider an intial value problem for a hyperbolic system of conservation 

laws 

Q. + F(Q) = 0 (2.1a) 

t X 

with 

Q(x,t=0) =■ Q q (x) (2.1b) 
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where the initial data Q 0 ( x ) are assumed to be piecewise smooth functions 

which are either periodic or of compact support. The state vector 

Q ■ (q^ ,^ 2 * • * ' >q m > ' * * >q M ) * s a c °l umn vector of M unknowns and the flux 

F(Q) is a vector-valued function of M components. The system in (2.1a) is 

hyperbolic in the sense that the MxM Jacobian A(Q) * 3F/3Q has M real 

eigenvalues X . The relation between the eigenvalues X and the Jacobian 
in m 

A(Q) is given by 


M 


A(Q) 


y r X l 
, m m m 
m*l 


( 2 . 2 ) 


where the columns r ffl are linearly independent right eigenvectors and the 
rows give an orthonormal set of left eigenvectors (i.e., Er^ i =* 6^). 

Let Qj and Qj+i denote two piecewise constant states of a Riemann 

problem. These two states are separated from each other by one or more wave 
families, which can be a shock, an expansion fan, or a contact discontinuity. 
Each wave family m(=l ,2 , . . . ,M) is propagated with its own wave speed which 
is given by the corresponding eigenvalue. In order to apply this wave model 
approach (cf. [16]) to equation (2.1a), consider (2.1a) for a computational 
cell with its centroid at Xj and its two interfaces at Xj ± y . The flux 
F^ i j is defined by 


F . i/ = F . - (Er X a ) . , i, m F . . - (ir X + a )., 
j+ Hi j v m m m'j+ Hi j+l v mm m'j+ Hi 


(2.3) 


where a . . i . = l .. i. A.,i. Q with A i.() - () ... - () . defines the 
m , J+ Hi m,J+V2 1+ Hi J+V2 J+l J 


strength of the m-th wave and X* * (lx |±X )/2. 

m 4 ® 4 m 
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Following Roe [8] , a mean value matrix 


A ( Q j ,Q j+l) ( A A )j + 1 / 2 ~ ^ Zr m X m A m~ Zr m X m ^nJj+V; 


+ ~7 T 7 


with the properties 


(i) A(Q,Q) - A(Q) 


(li) - 4j+1/2 F 


is utilized to write (2.3) as 


F.. i, - F, - Aj T/ A^Q - F . +1 - A+i / A j+V Q. 


*j ' 


(2.4) 


Using (2.4), a first-order accurate finite-difference approximation to (2.1) 
is easily constructed 


4Q ” - (t/h) ( V v 2 V v 2 Q ' a j- v 2 V v 2 Q)n ' 0 <2 • 5) 


where A() n = O^ 1 - () n , t is the time step size, and h - x ^ - x^_ y^ , 

An approximate solution computed with (2.5) tends to the exact solution for 
small data due to property (i). Property (ii) ensures that the approximate 
solution is exact if Qj and Qj +1 are connected by a single discontinuity, 
regardless of its size. This becomes evident by comparing (ii) with the jump 
condition 

v Vv 2 9 - VV 
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where V , 

Vv 2 9 


the discontinuity speed, is an eigenvalue of 
is the corresponding eigenvector. 


A< VW 


and 


Remark 2.1 : The derivation of the scheme in (2.5) tacitly implies that we 

have to deal with a Riemann problem at each and every interface. Then, the 

assumption of piecewise constant states Qj is equivalent to replacing the 

actual function Q(x) in x. i. < x < x i, with its cell average 

72 J+ >2 

i 

Q. * — J Q dx'. Consequently, the scheme in (2.5) and its higher-order 

J x j- v 2 

variants as discussed in this paper give approximations to cell averages of 


the exact solution rather than to its point values at the centroids x 


y 


Remark 2.2 : The scheme in (2.5) is not entropy satisfying [6] and admits 
stationary expansion shocks as its steady solution. Although this can be 
readily rectified [2,3], these entropy considerations will be ignored here for 
the sake of a simplified presentation. 


UNIFORMLY SECOND-ORDER ACCURATE TVD SCHEMES FOR SCALAR CONSERVATION LAWS 

For clarity we first consider the simplest possible case of a linear 
scalar equation 

q + aq^ “0, a>0. (3.1) 

The assumption a>0 does not impede generality; the contrary case is easily 
treated by symmetry. Applying the scheme in (2.5) to (3.1) yields 
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+ V W” " °> 


v * a(x/h) 


(3.2) 


where h is restricted to be constant, although results for irregular grids 
follow in a similar manner. A general second-order formulation of (3.2) reads 


Aq° + v[l+0.5(l-v)(4»"-<|»j_ 1 )]A :j _l /2 q n * 0 


(3.3) 


As usual [5,9,10,12,14], the quantities ^ are defined as functions of the 
ratio of two consecutive gradients 


, . n , . n N ±1 


for 4 j± v 2 1 * 0 


£or V'4’ ' 0 


(3.4) 


The scheme in (3.3) is conservative when the quantities <|>j satisfy Sweby's 
symmetry property [12] (i.e., <Jk * r^*). Christian the quantities 

<j>* with stencil selectors since they determine the underlying difference 
stencil of a particular finite-difference scheme. Choose 4 ^ = r^ and 

= const =*■ 1, and (3.3) recovers the second-order Lax-Wendroff scheme [7]. 

For * const « 1 and 4 >+ * r+, (3.3) is identical with the second-order 

Warming-Beam scheme [17]. Both linear schemes exhibit spurious oscillations 

at discontinuities of q(x) since they are not TVD [12] (i.e., not 

monotonicity preserving [2]). We turn now to find stencil selectors q>(r ^ ) 
which make (3.3) both TVD and uniformly second-order accurate. 


For a general scheme written in the form 
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(3.5) 


*£ - (c 


j j+ V2 A j+ V2 q C j" V2 fl j" V2 


- ct u. A, 1, q) n = 0 


where c*i. are data dependent (i.e., functions of the set {q 11 } = 
j + 72 n 


{,... ,q*j_^ .q^ , q ,•••}) , it is easily shown [2] that sufficient conditions 
for it to be TVD are the inequalities 


0<t rt and VV 2 + V V l - 


(3.6) 


A comparison of (3.3) and (3.5) yields 


c^. 1 / “ 0 and c j+ l/ 2 = v [ 1+0 * 5 ( 1_v ) < ♦ j+l"^ ] * (3 * 7) 


The scheme in (3.3) is TVD if 


-2/ ( 1-v) < - (f)^ < 2/v, 0<v<l 


(3.8) 


An infinite multitude of functions <j>^ = <j>(r*) satisfy condition (3.8). To 
derive one that is useful for the construction of uniformly second-order TVD 
schemes, additional constraints are required. Two of those constraints are 


<j>(0) = 0 
(f)(1) - 1- 


(3.9a) 

(3.9b) 


Condition (3.9a) is necessary to make (3.3) TVD [2], and condition (3.9b) 
immediately follows from Sweby's symmetry condition (i.e., <j>"t - rttj^), w ^ich 
makes scheme (3.3) conservative. 
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Having earlier established the notion of stencil selectors, makes it 
quite straightforward to define further bounds on <f>*. Second-order central 
differences do not introduce any dissipation and they have less dispersion 
than fully one-sided second-order differences. This suggests to control <}>* 
such that the scheme in (3.3) recovers the Lax-Wendroff scheme whenever there 
is no conflict with the TVD conditions in (3.8) and (3.9a). A function <j>(r^j) 
that satisfies these conditions is an antisymmetric extension of the min mod 
limiter 

<J>(r*) * max(Y,min(r*,l)), y - - 1 (3.10) 


which, for y m 0, is the conventional min mod limiter [2,4,5,9,12,18]. The 
conventional min mod limiter, like various other TVD limiters [9,12], returns 
<j>(r*<0) ■ 0 which reduces the scheme in (3.3) to a first-order upwind scheme 
even at smooth extrema of the solution q(x). Applying the antisymmetric min 
mod limiter in (3.10) to the scheme in (3.3) shows that second-order accuracy 
is maintained even across extrema (i.e., for r* < 0). Since the scheme 
becomes TVD of uniform accuracy, we call it a TVDUA scheme. It is interesting 
to note that equation (3.10) admits negative values for $(r*) ; a phenomenon 
incompatible with the notion of "antidiffusion" [9,12], which is often used to 
devise second-order TVD schemes . 


Remark 3.1 : In the remainder of this paper we will speak of first-order TVD 

schemes when 4>(r*)=const=0, second-order TVD schemes when <j>(r*) is defined 
by a conventional min mod limiter, and of TVDUA schemes when <£(r^) is 

defined by an extension of the conventional TVD limiter as in (3.10). 
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The extension of the TVDUA scheme to a nonlinear conservation law like 


q t + a(q)q x - o 


(3.11) 


is straightforward. Applying the scheme in (2.5) to (3.11), we obtain 


S ' ( VV 2 Vv 2 ’ • HVV 1 " ■ 0 


(3.12) 


where v^_]^ = a j- ly t/h is now a local CFL number with a" = (|a|±a)/2. 
Using the general form in (3.5), a conservative second-order extension to 
(3.12) which accounts for the local flow direction has the coefficients 

c 5+v 2 " VV 2 [ 1 " 0 - 5(1 'VV 2 , *j+i ] + °' 5 v j-V 2 (1 'VV 2 ) *j <3 - 13a) 

c j+ V 2 * VV 2 + °- 5 V3/2< 1 'V3/2>*j + i '(3.13b) 

The coefficients in (3.13) preserve conservancy whenever the limiters <J>* 
satisfy Sweby^s symmetry property [12] (i.e., (j>^ » and fc key satisfy the 

TVD conditions (3.6) when 


■ 2/ V V 2 ‘ Vi - +3 ‘ 2/ v v 2 


(3.14a) 


■ 2/(1 ‘VV 2 ) ‘ Vi - *j ‘ 2/( i'V 1 / 2 > 


(3.14b) 


These restraints of are satisfied for all 0 < v* i, <1 when <j> + is 

1 J+ 72 

chosen is already defined in equation (3.10). 
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UNIFORMLY SECOND-ORDER ACCURATE HIGH RESOLUTION SCHEMES FOR HYPERBOLIC SYSTEMS 

The TVDUA schemes are formally extended to hyperbolic systems in the 
usual field by field decomposition [1,2,4,5,9,12,18]. Defining characteristic 
variables: 

w * l m Q 

m m 

and assuming the constant coefficient case (i.e., A is a constant MxM 
matrix), equation (2.1a) is written as a set of M decoupled linear scalar 
conservation laws 

( w m ) *. + * ( w m ) Y * 0 (4.1) 

m t m m x 

A first-order upwind approximation to (4.1) reads 



l 


(T/h)(x; i}tl/ W n - X*l J . 1/ W m ) n - 0 


(4.2) 


The uniformly second-order accurate TVDUA scheme for linear scalar 
conservation laws as defined in (3.3) and (3.10) is straightforwardly 
applicable to (4.2): 


Aw n - {v~[l-0.5(l-v )(<T ...V .)]A.. i. w 
j 1 m L m m, j+1 r m,j j+ l h m 


- v + [ 1+0 . 5 ( 1- v + ) ( . .)]A. i, w } n = 0 (4.3) 

m L m m,j T m,j-1 J j - V? v 7 


with v* = X* x/h and 
m m 
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(4.4) 


*< r m, j> " "“H* min(r m, j* 1 ^ 


where 


( V'4VVv ! , : ) 


n*±l 


for A . . i , w n * 0 


®»J 


j±V 2 « 


for 4 j*V 2 w »' 0 


(4.5) 


The nonlinear system case is recovered in two steps [5]: first, the TVDUA 

scheme in (4.3) is modified as 


Aw 11 . - {v i, [l-0.5(l-v _ i,)(j) . .] 

®>J 1 ®»j+V2 L m,j+V2 ®,J+ 1 J 


+ 0.5 v . i, (1-v . i. )$ ,} n A . , i , w 11 

m,j-V2 m.r /2 “»j J+V2 ™ 


m 


+ {v + • 1/ [l*0.5(l-v + . i. H + . ,1 
1 ®»J-V2 L ®»J-V2 


+ 0.5v + ^ i, (l-v + i, )<f .} n A. i. w n - 0 
m, j+ V2 ®, 1+ V2 ®,J J J- ¥2 ® 


(4.6) 


to accommodate the nonlinear scalar conservation law case. The stencil 

selectors d>(r^ .) remain the same as in (4.4). Following a multiplication 
ni, J 

of (4.6) by the set of the right eigenvectors from the left, a 

conservative, uniformly second-order accurate high resolution scheme for the 
nonlinear systems case is easily derived 
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49 J 'l f r »,j + V 2 v »,j+ 1 /2t 1 '°- 5<1 ' v « > j + V 2 ) Vj + il 


+ 0.5 r . I, v .i, (1-v . i, )<t> .} n A aj_1/ w 11 

m,J-V2 m,J-V2 m,j-V2 m,J J 1+ L /2 m 

+ 7 {r , i, v + . i, [l-0.5(l-v + . 1 . )<J> + . . ] 

^ l m,j-V2 ®* J~ V 2 L m > j~ <2 m »J -1 

+ 0.5r i. v + 1 , (l-v + i, )<|> . } n A . h w n = 0 (4.7) 

m> J+ t 2 m»J+V2 m,J+V2 m,j J J- 72 * 


Although it is no longer necessarily TVD, the scheme in (4.7) gives 
highly resolved solutions in regions of smooth Q(x,t) while spurious 
oscillations are suppressed in regions of rapid changes in gradient. 


NUMERICAL EXAMPLES 

Consider the constant coefficient case 

(5.1a) 
(5.1b) 


q + q 

^t x 


with 


q(x,t=0) = sinirx, 0 < x < 2 


and periodic boundary conditions for an initial numerical test. The exact 
solution to (5.1) for any given time is easily computed by solving the 
characteristic equation 

q(x,t) ■ sin ir(x-t) (5.2) 
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Let the interval [0,2] be divided into (J-l) equidistant computational cells 
with their centroids being defined at 

x - 2( j-l)/ (J-l) , 1<J<J (5.3) 

Figure 1 shows comparisons of the exact solution (solid line) with three sets 
of numerical results, at t*2 with x/h a 0.8 and J=21. As expected, the 
first-order accurate TVD scheme (fig. la) approximates the exact solution 
rather poorly: a consequence of its hefty inherent numerical dissipation. A 
comparison of figures lb and lc demonstrates how the degeneracy to first-order 
accuracy at local extrema in the second-order accurate TVD scheme adversely 
affects the global accuracy in the neighborhood of extrema. This results is 
substantiated by the grid refinement study in figure 2 for J=21, 41, and 
81. The computations with the first-order TVD scheme converge with first- 
order accuracy in both the Lj- and the error norm. The TVDUA scheme is 

second-order accurate in Li as well as in L ; in between lies the second- 
order TVD scheme which is second-order accurate in L^ but only first-order 
accurate in L . 

oo 

Now we turn to examine numerical approximations to a nonlinear scalar 
problem defined by 

q t + (q 2 /2) x - 0 (5.4a) 

q(x,t=0) = a + 6 sinir(x+Y), 0<x<2 (5.4b) 

q(x=0,») = q(x=2,«) (5.4c) 
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For 0<t<t Q , the solution to (5.4) is smooth and is "exactly” computed by 
using, for instance, Newton-Raphson iterations to solve the characteristic 
relation 

q(x,t) ■ a+ 6 sinir(x+Y-qt) , 0<t<t Q (5.5) 

When t=*t Q 3 shock develops at x Q and moves with its shock speed q Q for 

some time before it starts interacting with the expansion wave which brings 

about a rapid decay of the solution. The values q Q> t Q , and x Q are 

2 2 

computed by requiring 3q/3x * » and 3 q/3x =0. We find that a shock 

develops at t = 1/Bir, x » 1— y+a/Bir> and that it moves for some time with 
o o 

its speed q (x ,t ) = a. 
o o o 

Figures 3, 4, and 5 present computations of (5.4) for a=l, 8-0.5, and 
Y*0 (i.e., q(x,0) = 1 + 0.5 sin irx); thus, a shock forms at t=2/ir. The 

calculations are carried out on meshes as defined in (5.3). In figure 3, the 
"exact" solution ("solid line) is compared with three calculations of 
increasing accuracy at t*0.3 with x/h-0.6 and J-21. Particularly the 

approximation of the "exact” solution in the region 1.5<x<2 reveals the 
superior resolution capability of the TVDUA scheme when compared to the first- 
and second-order accurate TVD schemes. The corresponding grid refinement 
study in figure 4 shows the convergence of the L, and the L error norms 
for J-21, 41, and 81. The solutions with the first-order TVD scheme converge 
with first-order accuracy in both the and the error norm. The 

second-order TVD scheme proves to be second-order accurate in the L-^ norm, but 
it is barely more accurate than the first-order TVD scheme in the L error 

oo 

norm. The TVDUA scheme is clearly second-order accurate in the norm and 
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almost second-order accurate in the L norm. The effect of error 

00 

accumulation in time becomes apparent when the initial differences between 
"exact" and numerical solutions for t»0.3 in figure 3 are compared with the 
corresponding results in figure 5 for t«2/ir (i.e., when the shock develops). 
The first-order accurate results deviate considerably from the "exact" 
solution. The second-order TVD scheme captures the forming shock quite 

decently while the TVDUA scheme "nails" the "exact" solution. 

Finally, the TVDUA schemes are put to test when applied to nonlinear 
hyperbolic systems. It is difficult to find a suitable test case since the 
extra resolution power of our TVDUA schemes like that of the UNO schemes by 
Harten and Osher [4] manifests itself best across smooth extrema of a 
solution. Consequently, we found the commonly chosen shock tube problems of 
Lax or Sod [2,3,4] not very illustrative because their solutions are just 
constant states separated by waves. In order to at least achieve a smooth 
variation of states, we consider a one-dimensional flow in a duct of variable 
cross-sectional area B(x) which is described by 

Q + F - H 
x t x 

T 

Q * (P.pu.E) 

F - (pu,(pu 2 +p), (E+p)u) T 

H - (0,p(dinB/dx),0) T (5.6) 

with the density p, the velocity u, the pressure p, and the internal energy 

2 

per unit volume E = pe+pu /2 (e: specific internal energy). The eigenvalues 

of the matrix A = 9F/9Q are 
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X ^ * u-c, * u, X^ ■ u+c (c: speed of sound). 

The corresponding sets of orthonormal left and right eigenvectors are given 
elsewhere (cf. [2,8,18]). Using Roe's averaging [8] to compute the mean value 
matrix ^(Qj'Qj+l)* ttie a PP roximate Riemann solver as described in Section 2 
is employed to obtain a first-order upwind approximation to (5.6). The second- 
order versions of that baseline scheme are constructed by formally extending 
the TVD and TVDUA concept to hyperbolic systems as described in Section 4. 

The geometry of the divergent duct which is depicted in figure 6, is 
taken from a rather common benchmark test [2,10,11,18]. Specifying 

p, - 0.502, u = 1.299, e = 1.897 
in in in 

at the inflow cross section, a stationary compression shock forms at x = 4.816 
[13,181 when p = 0.776 at the outflow cross section. The shock causes a 
problem in quantitatively determining the accuracy of the employed schemes, 
since the notion of accuracy being strongly coupled with the notion of 
truncation error in finite-difference methods becomes immaterial for 
nondifferentiable solutions. Despite this shortcoming, the results in figures 
7 and 8 demonstrate the superior resolution power of the TVDUA schemes when 
compared to the first- and second-order TVD schemes. All computations are 
carried out on equidistant meshes with just 21 grid points and with T/h=0.2. 
The "exact” solution (solid line) and their numerical approximations are 
plotted in terms of density (fig. 7) and velocity (fig. 8) as a function of x. 
The first-order TVD scheme does not only smear the shock, but it also 
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noticeably misses the "exact” solution for 2.5<x<7.5: a range consisting of 

ten (!) grid points. The approximate solution computed with the second-order 
TVD scheme feels the ripples of the shock in a narrower regime with 3<x<6 
which translates to 7 grid points. Both TVD schemes resolve the shock with 
two intermediate zones, whereas, the TVDUA scheme captures the shock in just 
one interval. The TVDUA results also follow the "exact" solution much closer 
ahead and aft of the shock. 


CONCLUDING REMARKS 

A class of new explicit, uniformly second-order accurate high resolution 
schemes, which are TVD in the scalar, or linear systems, case, has been 
constructed, analyzed, and verified in several test calculations. Compared to 
conventional high resolution schemes which switch back to first-order accuracy 
at extrema of the solution, the novel schemes here called TVDUA schemes, yield 
a gain in accuracy at virtually no computational extra cost. The gain in 
accuracy becomes most apparent when solutions with smooth extrema are 
approximated as in the test calculations presented for scalar conservation 
laws. When dealing with essentially discontinuous solutions, the TVDUA 
schemes produce nonoscillatory solutions, they resolve shocks with just one 
interval, and they sustain their extra resolution capability almost unaffected 
even in the the neighborhood of a strong singularity. Future work will 
concentrate on the extension of the novel TVDUA schemes to multi-dimensional 
problems. 
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Figure 2. Error norms for numerical solutions to an one-dimensional, 


linear wave equation: t=2, T/h=*0.8. 
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Figure 5(b) 

Figure 5. Exact and numerical solutions for an one-dimensional, wave 
equation problem at the onset of a shock formation: t~2/n, 

x/h-0.6, 21 grid points. 
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